{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 367,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from numpy import sqrt,mean,log,diff\n",
    "import QuantLib as ql\n",
    "from pandas_datareader.data import Options\n",
    "import pandas_datareader.data as web\n",
    "import datetime "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 368,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "opt = Options('spy', 'yahoo')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 369,
   "metadata": {},
   "outputs": [],
   "source": [
    "expiration_dates = [ql.Date(i.day, i.month, i.year) for i in opt.expiry_dates]\n",
    "expiry_index = 14 # choose the contracts expire in 4 months\n",
    "data = opt.get_call_data(expiry=opt.expiry_dates[expiry_index])\n",
    "strikes = list(data.index.get_level_values('Strike'))\n",
    "premium = list(data['Last'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 534,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "datetime.date(2017, 11, 17)"
      ]
     },
     "execution_count": 534,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt.expiry_dates[expiry_index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 370,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "day_count = ql.Actual365Fixed()\n",
    "calendar = ql.UnitedStates()\n",
    "calculation_date = ql.Date(opt._quote_time.day,opt._quote_time.month,opt._quote_time.year)\n",
    "spot = opt.underlying_price\n",
    "ql.Settings.instance().evaluationDate = calculation_date\n",
    "dividend_yield = ql.QuoteHandle(ql.SimpleQuote(0.0))\n",
    "risk_free_rate = 0.01\n",
    "dividend_rate = 0.0\n",
    "flat_ts = ql.YieldTermStructureHandle(\n",
    "    ql.FlatForward(calculation_date, risk_free_rate, day_count))\n",
    "dividend_ts = ql.YieldTermStructureHandle(\n",
    "    ql.FlatForward(calculation_date, dividend_rate, day_count))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 518,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# dummy parameters\n",
    "initial_var = 0.2; rate_reversion = 0.5; long_term_var = 0.2; corr = -0.75; vol_of_vol = 0.2;\n",
    "# initial_var = 0.2; rate_reversion = 0.15; long_term_var = 0.6; corr = -0.75; vol_of_vol = 0.2;\n",
    "process = ql.HestonProcess(flat_ts, dividend_ts, \n",
    "                           ql.QuoteHandle(ql.SimpleQuote(spot)), \n",
    "                           initial_var, rate_reversion, long_term_var, vol_of_vol, corr)\n",
    "model = ql.HestonModel(process)\n",
    "engine = ql.AnalyticHestonEngine(model) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 519,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "heston_helpers = []\n",
    "date = expiration_dates[expiry_index]\n",
    "for j, s in enumerate(strikes):\n",
    "    t = (date - calculation_date)\n",
    "    p = ql.Period(t, ql.Days)\n",
    "    sigma = premium[j]\n",
    "    helper = ql.HestonModelHelper(p, calendar, spot, s, \n",
    "                                  ql.QuoteHandle(ql.SimpleQuote(sigma)),\n",
    "                                  flat_ts, \n",
    "                                  dividend_ts)\n",
    "    helper.setPricingEngine(engine)\n",
    "    heston_helpers.append(helper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 520,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lm = ql.LevenbergMarquardt(1e-8, 1e-8, 1e-8)\n",
    "model.calibrate(heston_helpers, lm, \n",
    "                 ql.EndCriteria(500, 50, 1.0e-8,1.0e-8, 1.0e-8))\n",
    "long_term_var, rate_reversion, vol_of_vol, corr, initial_var = model.params()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 521,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "long_term_var = 0.191762, rate_reversion = 0.000001, vol_of_vol = 0.215442, corr = -0.817388, initial_var = 0.198778\n"
     ]
    }
   ],
   "source": [
    "print \"long_term_var = %f, rate_reversion = %f, vol_of_vol = %f, corr = %f, initial_var = %f\" % (long_term_var, rate_reversion, vol_of_vol, corr, initial_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 522,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mc_heston(option_type,S0,K,T,initial_var,long_term_var,rate_reversion,vol_of_vol,corr,r,num_reps,steps):\n",
    "    \"\"\"\n",
    "    option_type:    'p' put option 'c' call option\n",
    "    S0:              the spot price of underlying stock\n",
    "    K:              the strike price\n",
    "    T:              the maturity of options\n",
    "    initial_var:    the initial value of variance\n",
    "    long_term_var:  the long term average of price variance\n",
    "    rate_reversion: the mean reversion rate for the variance\n",
    "    vol_of_vol:     the volatility of volatility(the variance of the variance of stock price)\n",
    "    corr:           the correlation between the standard normal random variables W1 and W2\n",
    "    r:              the risk free rate\n",
    "    reps:           the number of repeat for monte carlo simulation\n",
    "    steps:          the number of steps in each simulation\n",
    "    \"\"\"\n",
    "    delta_t = T/float(steps)\n",
    "    payoff = 0\n",
    "    for i in range(num_reps):\n",
    "        vt = initial_var\n",
    "        log_st = log(S0)\n",
    "        for j in range(steps):\n",
    "            w1 = np.random.normal(0, 1)\n",
    "            w2 = corr*w1+sqrt(1-corr**2)*np.random.normal(0, 1)\n",
    "            vt = (sqrt(vt) + 0.5 * vol_of_vol * sqrt(delta_t) * w1)**2  \\\n",
    "                 - rate_reversion * (vt - long_term_var) * delta_t \\\n",
    "                 - 0.25 * vol_of_vol**2 * delta_t\n",
    "            if vt < 0: vt = 0.00\n",
    "            log_st = log_st + (r - 0.5*vt)*delta_t + sqrt(vt)*sqrt(delta_t)*w2\n",
    "        st = e**(log_st)\n",
    "        if option_type == 'c':\n",
    "                payoff += max(st - K, 0)\n",
    "        elif option_type == 'p':\n",
    "                payoff += max(K - st, 0)\n",
    "        \n",
    "    return (payoff/float(num_reps)) * (exp(-r*T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 523,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.383434958491186"
      ]
     },
     "execution_count": 523,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mc_heston('c',spot,strikes[20],t/365.0,initial_var,long_term_var,rate_reversion,vol_of_vol,corr,0.01,100,1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 524,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36.119999999999997"
      ]
     },
     "execution_count": 524,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "premium[20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 525,
   "metadata": {},
   "outputs": [],
   "source": [
    "heston = []\n",
    "\n",
    "for i in range(len(strikes)):\n",
    "    heston.append(mc_heston('c',spot,strikes[i],t/365.0,initial_var,long_term_var,rate_reversion,vol_of_vol,corr,0.01,100,100))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 529,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = [(heston[i]-premium[i])**2 for i in range(len(premium))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 530,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8885.1013305184661"
      ]
     },
     "execution_count": 530,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(diff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 531,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(110.46562812345368, 127.83),\n",
       " (110.6988943943971, 112.12),\n",
       " (55.926776121382609, 55.710000000000001),\n",
       " (52.981820291517408, 54.030000000000001),\n",
       " (54.164283840319435, 53.340000000000003),\n",
       " (54.991136360501223, 52.719999999999999),\n",
       " (41.216283956629468, 51.090000000000003),\n",
       " (47.687719919883207, 50.340000000000003),\n",
       " (44.545781271634148, 49.32),\n",
       " (47.563439672575591, 48.340000000000003),\n",
       " (44.635711823762009, 47.340000000000003),\n",
       " (46.252466035738571, 46.170000000000002),\n",
       " (47.574492181620364, 45.18),\n",
       " (41.653315658070859, 43.310000000000002),\n",
       " (47.604806657885561, 43.229999999999997),\n",
       " (40.894480554924115, 42.799999999999997),\n",
       " (38.539091518531059, 42.25),\n",
       " (43.303139080062273, 40.829999999999998),\n",
       " (36.234352000924453, 39.880000000000003),\n",
       " (39.257917188012463, 38.469999999999999),\n",
       " (44.853854689174561, 36.119999999999997),\n",
       " (39.640305408524711, 36.310000000000002),\n",
       " (34.129177338407089, 36.359999999999999),\n",
       " (36.438197093590269, 35.5),\n",
       " (39.905750087595855, 34.130000000000003),\n",
       " (32.778166732302317, 32.649999999999999),\n",
       " (43.333301445006789, 31.969999999999999),\n",
       " (36.191610071798877, 30.550000000000001),\n",
       " (32.439741216413196, 31.300000000000001),\n",
       " (26.253346286101117, 29.34),\n",
       " (40.260711758568831, 28.760000000000002),\n",
       " (30.35683234367302, 27.449999999999999),\n",
       " (26.437940074957758, 27.579999999999998),\n",
       " (28.483162465814967, 25.710000000000001),\n",
       " (25.988006256896163, 25.809999999999999),\n",
       " (27.515272593610401, 23.98),\n",
       " (25.183307114602915, 22.760000000000002),\n",
       " (23.158436475438911, 21.850000000000001),\n",
       " (33.342453057590774, 21.100000000000001),\n",
       " (27.131945603771868, 20.82),\n",
       " (23.601047412852459, 18.010000000000002),\n",
       " (25.760858340591355, 17.719999999999999),\n",
       " (23.691135308980339, 16.440000000000001),\n",
       " (22.185597107961662, 15.16),\n",
       " (20.277352259807657, 14.58),\n",
       " (27.451461967610275, 13.720000000000001),\n",
       " (31.744386996748887, 12.76),\n",
       " (26.49106582685053, 11.970000000000001),\n",
       " (21.003755215426228, 11.25),\n",
       " (19.35608163937891, 10.279999999999999),\n",
       " (17.351500422342902, 9.5999999999999996),\n",
       " (22.659139445207355, 8.7899999999999991),\n",
       " (23.409228489912628, 8.0700000000000003),\n",
       " (20.961754171171588, 7.5899999999999999),\n",
       " (16.075588910706614, 6.6600000000000001),\n",
       " (20.825896797061727, 6.0),\n",
       " (18.675398030756703, 5.5),\n",
       " (19.026462754554451, 4.79),\n",
       " (19.417525830223067, 4.21),\n",
       " (20.407591273285487, 3.7000000000000002),\n",
       " (17.174337684195297, 3.1800000000000002),\n",
       " (17.391313889375883, 2.71),\n",
       " (16.050599868535851, 2.2000000000000002),\n",
       " (14.09485225440975, 1.8899999999999999),\n",
       " (18.208532190234333, 1.48),\n",
       " (17.350492587503762, 1.1699999999999999),\n",
       " (15.66130115227647, 0.96999999999999997),\n",
       " (18.310204566097372, 0.71999999999999997),\n",
       " (15.244762462507955, 0.55000000000000004),\n",
       " (12.05896923742357, 0.41999999999999998),\n",
       " (16.549370172554656, 0.32000000000000001),\n",
       " (15.127520335730379, 0.25),\n",
       " (10.679068848118462, 0.19),\n",
       " (17.091513222601169, 0.16),\n",
       " (15.164841725875474, 0.13),\n",
       " (12.162076252506601, 0.12),\n",
       " (11.179740065544316, 0.10000000000000001),\n",
       " (15.126628346299601, 0.089999999999999997),\n",
       " (10.827121868063745, 0.080000000000000002),\n",
       " (7.955281416676554, 0.080000000000000002),\n",
       " (11.584285727318449, 0.070000000000000007),\n",
       " (6.8753197768030745, 0.070000000000000007),\n",
       " (10.195241428465291, 0.070000000000000007),\n",
       " (11.098792968645823, 0.080000000000000002),\n",
       " (10.97397900549743, 0.059999999999999998),\n",
       " (8.8726648972088906, 0.059999999999999998),\n",
       " (11.791012392977612, 0.040000000000000001)]"
      ]
     },
     "execution_count": 531,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zip(heston, premium)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
